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ABSTRACT 


A  practical  method  for  computing  the  conditional  expectation  of  a 
polynomial  in  the  components  of  a  multivariate  normal  random  variable 
X,  when  X  is  restricted  to  a  subset  of  »?,  is  given.  This  method 
makes  the  application  of  certain  missing  data  techniques  possible  in 
cases  there  repeated  numerical  integration  is  not  feasible. 
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APPROXIMATING  CONDITIONAL  MOMENTS  OF  THE  MULTIVARIATE 
NORMAL  DISTRIBUTION 


By 

Joseph  G.  Deken 


1.  Introduction. 

The  conditional  moments  such  as  E (X^ )  of  a  multivariate  normal  random 

u 

variable  X  =  (X^Xg,. .  •  ,X^),  when  X  is  restricted  to  a  subset 
Ad**  are  not  readily  obtained  numerically,  since  the  required 
integration  in  p -dimensions  is  time-consuming  except  for  very  small 
p.  These  conditional  moments  are  of  interest,  for  example  in  the 
derivation  of  E-M  estimates  in  missing  data  problems  (Dempster,  Laird 
and  Rubin,  1977).  We  present  here  an  efficient  approximation  scheme  for  these 
moments,  which  makes  the  computation  practical  for  moderately  large  p. 


For  convenience  of  description,  we  restrict  attention  to  sets 

A  of  the  form  I.  x  i  x  • • •  x  I  ,  where  all  the  I .  are  intervals, 
1  2  p'  j 

but  the  approach  is  more  general,  as  indicated  below.  We  start  by 
observing  that  the  approximations  to  the  ratio 


i-tJ 


2<t2 


dx 


B+tr 

1-tJ 


2 <r 


dx 


obtained  by  the  first  terms  in  a  Taylor  series  around  t  =  0  may  be 
written  as  polynomials  in  p: 
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The  conditional  expectation  of  a  polynomial  P(X)  =  a^+a^X+a^X2 +•  •  • 
in  the  normal  random  variable  X  is  thus  approximated  by  replacing  P 
by  a  polynomial  Q(p)  in  the  mean  of  X,  where  the  transformation 
T(s,t,cr  ):  P  Q  is  linear.  Since  the  transformation  V:  c  ■*  b  defined 

by 


is  also  linear,  and  the  conditional  mean  of  Xp  given  X^, . . .  ,X^  is 
of  the  form  (g-HiX^^)  where  g  is  linear  in  X1,...,Xp  2,  any  approxi¬ 
mation  to  EX^  which  is  a  polynomial  in  the  mean  of  Xp  will  produce  a 
polynomial  in  (X^ . . .  ,X^_^)  when  Xp  is  conditional  on  XL,...,Xp_1. 

The  conditional  expectation  in  3R*5  may  thus  be  accomplished  in  2p 
steps,  by  applying  T  and  then  V  in  succession  p  times,  to  obtain 

E(',Xl"*XP-l)  WE('IXl---Xp-2)  "  *“  "  E(*  |X1}  "Jfc  Wl  * 

Computationally,  this  process  requires  2p  matrix  multiplications  and  is 
thus  practical  for  moderate  p. 


2.  An  Example. 

The  following  simple  example  will  serve  to  establish  some  Ideas 
and  notation.  Let 

(u-x)2 

2  S+tf  V  "  c*2 

L  (o-  ,p,s,t)  =  x  e  31  dx  . 

s-t 

We  are  concerned  with  the  conditional  expectation  I^/Iq  of  X*,  where 
X  is  a  normal  random  variable.  The  sum  of  the  first  terms  of  a  T^lor 
series  for  this  ratio  about  t  =  0  is  of  the  form 


(1) 


2 


a=o 


*  /  2  \ ,  2CX 

CltfX  s  * 


but  also  may  be  represented  by  identifying  the  coefficients  of  4“  in  (l) 
as 

2£  (  2  *  \  a 
L  cfv'0’  • 


a=dO 


For  example,  (letting  v  =  or  ) 


h/To  • 


■  (n-a)t2  ,  gv(a-ii)  *  (s-u)5 
5v  1.5V5 


•  "(1 


2  2  4  2  k  2  2k  ?  li 

-lstv+2tv+st  X  ,  ,15t  v  -2t  v-3s  t\ 

57  0  ( — 57 - )u 

U  U 

/St  x  2  /  t _ \  it 

v - 5-;u  -  ( - r)u  +,‘' 

15v 
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j.  Multiple  Integration. 


If  (X1,X2 )  are  bivariate  normal,  the  conditional  distribution  of 

2 

X2,  given  X1,  is  TI(m2-Hd21(X1-m1),  cr2 ^3.  The  approximation  procedure 


given  above  yields  (writing  p  for  the  vector  of  coefficients  of  P,  and  with  C 

**  k 

the  matrix  whose  j,ktb  element  is  the  coefficient  of  4  in  the  approximation 


of 


rA 


): 


E(P(X2)|X2e(s2-t2,s2+t2),X1)  =  (^C(s2,t2,a2>1))(42+b21(X1-Ml)) 


This  last  is  of  course  a  polynomial  in  X^,  with  coefficients 


:=  b21  ^  ^£C^S2,t2,<T2.1^P^a^>i2"b21(il^  ' 

p  >  ut 


so  that 


E(p(X2)|X2e(s2-t2,s2H2),  XLe(s1-t1,s1+t1))  = 

(qC  (s^,  t^,o"i ) )  (^i )  . 

Other  cases  (higher  dimensions,  multivariate  polynomials)  are  treated 
similarly,  but  it  will  be  best  to  introduce  some  notation  at  this  point. 
If  p  is  a  vector  of  coefficients,  then 

SUBS(a,b,g)  is  a  vector  of  coefficients, 

(SUBS(a,b,p) ) .  :=bJ  £  a*"J  • 

0  k  >  j  0 

(i.e.  SUBS(a,b,p)  is  the  vector  of  coefficients  of  y^  when 
a+by  is  substituted  for  x  in  p(x)). 

U 


With  the  aid  of  a  symbolic  computation  system  such  as  MIT' s  MACSYMA, 
the  coefficients  of  the  transformation 


a 


a 


; 


as  functions  (v,  s,t)  may  be  computed  for  much  larger  values  of  k 

and  a  then  would  be  practical  by  hand.  (A  partial  table  of  values  c 
is  given  in  ap  endix  A  of  this  paper.  Card  versions  consisting  of  FORTRAN 
assignment  statements  "  c  (l,j)  =  •*•  "  may  be  obtained  on  request  from 
the  author. )  As  subsequent  examples  show,  the  approximation  of  in 
the  box  1^  x  x  •  •  •  I  involves  cffla  for  higher  values  of  m  than 
k,  depending  on  the  order  of  the  underlying  Taylor  series.  Fortunately, 
the  approximation  based  on  only  a  few  terms  is  quite  accurate,  as  evidenced 
by  the  following  example:  Since  the  integral 


.  (s-t)2 

xe  2  dx  =e  2 


(s+t)2 

2 


the  Taylor  series 


Ij. (l,o,s, t)  t2  Sts2  k. 

^  (1,0,.,  t)  -  S(1  T'TT*  ) 


gives 


(s-t)  _  (s+t)* 


I  (1,0, s,t)  =  <T>(s+t)-$(s-t)  » 


2  2 

e  -  e _ 

t2  ,  2+s2  .4 ■ 


•  (l  -  y  +  *  ) 

The  approximate  values  obtained  for  *(x)  =  .5  +  $(x)  -  *(0)  by  this 


method  agree  with  tabled  values  for  four  decimal  ’ilaces  for  0  <  x  <  1.1. 


r 


INT(s,t,v,p)  is  a  map  from  the  scalars  s,t,v  and  the  vector  p  to  the 
vector  pC(s,t,v),  i.e.  matrix  postmultiplication,  where  we  take 

m 

L /ln  ~  £  c  as  defining  C. 

K  0  a=0  J 


We  will  also  consider  INT  as  a  function,  so  that  if  INT  =  (aQ,a1,. . . ,am), 

then  (INT ) (x )  =  aQ+a^X  +•  •  •  +  a^X0  • 

With  the  above  notation,  we  can  define  multiple  integrals  in  a  fairly 

compact  fashion.  Hie  multivariate  normal  distribution  may  be  parameterized 

by  a  vector  of  means  *  ♦  >  Up)*  the  conditional  variances  (V^,...,Vp) 

defined  by  VI  =  Var(X.  ),  V.  =  Var(X . [X., ,. . .  ,X.  ^  ),j  >  1,  and  the  regression 

coefficients  b  7 ,b  0,b,1 , . . . ,b  b.  defined  by 

<z±  jL  pjP“J-  P-*- 


E(X 


X, , . . 


J-l 

•^J.l)  =  £ 

J  x  i=l 


en  Vi 


One-dimensional  integration  is  done  directly: 


(P»s,t,n,v)  4-  (lNT(s,  t,v,p) )  (u)  . 

Further  integrals  (12,13,...)  in  higher  dimensions  are  defined  recursively, 
e.g.  in  two  dimensions: 

I2(s1,t1,s2,t2,n1,M2,v1,v2,b21,p)  :  = 

(INT  (s^^^, SUBS  )M2_b21^1'b2l,  ^(Sg,^,  v2,£)  ) ) )  * 

and  from  this  we  obtain  the  two-dimensional  integral  as 

^  ^2*  ^1*  V1,V2,^21,£^  ^  ^^1  ^  " 
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In  general,  IK  =  INT(. . .  SUBS (. . . I(K-l) ) ).  We  have  restricted  attention 


here  to  sets  of  the  form  I.  x  I  x  ...  I  where  the  I.  are  intervals. 

1  2  p  j 

For  more  general  sets  A,  a  similar  scheme  would  work  provided  the 
section  of  X.  e  given  X1  •••  ^  is  an  interval  with  endpoints 

idlieh  are  polynomials  in  •  •  •,  X^_1< 
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k.  A  Numeric  Implementation. 


The  following  example  shows  how  the  three-dimensional  conditional 

expectation  approximation  may  be  defined,  using  only  vector  and  matrix 

arithmetic  with  numeric  arguments.  The  polynomial  to  be  approximated 

is  presumed  to  be  of  degree  at  most  two,  and  perhaps  bivariate  in 

and  X^.  (This  case  is  general  enough  for  all  the  means,  squares,  and 

products  EX^EXg,...  EX^,EX^, ...  EX^X^ ...  EXgX^).  In  fact,  we  describe 

here  the  approximate  values  of  E  (X3),  E(x|),  E (X  2Xj),  since  other 

expectations  reduce,  by  interchange  of  variables,  to  these  three.  The 

4 

approximation  used  will  include  powers  of  t  up  to  t  inclusive,  i.e. 


»  £  c  (s,t,v)“  . 


^  6 

The  functions  (c.  }"  may  be  found  in  appendix  A). 

ja  a=o 

To  approximate  E(X^),  the  integration  on  X^  is  first  carried  out, 


yielding 


£  clot(  Wv3  ^M3"Mlb5l"M2b32+Xlb51+X2b32^ 


3  a  a-£  a  _ 

=  a^  cb(83*Vt3^  ^  (b31Xl ''  (b32x2)  (,a3"^lb3l'M2 

3  3  -o:  3  7 

a5)  P5)  7  J^(°t^)Cl>(S3,t3'V3)(M3“Mlb3l"M2b32)7  (  ^)b31b32  * 


3  3-Q  a  p 

The  above  expression  is  of  the  form  £  £  r^,  X.X  ,  where 

a=0  P=0  ^  1  d 


raB  "  b31b32  *  ^3">1lb3l‘M2b32)7  ^  ^ 


Integration  on  X^  then  yields: 


3  3-a  „  3 

E  £  rofl  (  £  ^VVV 
a=o  p=o  ^  7=0  P7 


(M2~^lb21+b21Xl ^  ^ 


6  a  &A3  3^6 


=  E  sa^i  >  ^ere  s  =  £  £  r  £ 

—  a  6=(a-2)V0f&)  8p7=a-6^-6  21 


a=o 


C^7  ^S2,t2,V2^  ’  ^(i2“Mlb21^ 


7- (a- 6) 


and  the  final  form  is 


a£>  sa  £  *<VVV*S 


The  approximation  of  X^  is  carried  out  in  the  same  fashion: 

^  ‘i  s“  X  - 


a  A3 


3-B 


Si  '  M?.5)»o  e5>  r»  rJ.6(°-6)bS%r<a2't5’T2)  ‘ 


/  -u  %7-(0-B) 
^_Mlb2l'  } 


rag  ~  b31b32  7 ^c27^s3,t3,V3^^3"Mlb3l"M2b32^ 


(Note  that  the  only  charge  is  the  substitution  of  )  for 

cl?(s3’t3’v5,)- 
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To  approximate  XgX, , 


the  first  approximate  integration  is  exactly 


■4  1 


as  in  approximating  Xy  but  the  resulting  sum  is 

X2^  E  E  ra(3XlX2^  ’ 

a=o  p=o  H 

i.e.  we  have  r^  =  a  second  integration  gives 


3  W*  «  3 


JE^  rapXl^^  %y  ^S2^t2,V2^^2-b21Ul+b21Xl^  ^ 


so  the  final  result  is 


where 


c£)  S“  p?0  W^VV^l  ’ 


a  A3  4-8 


a  5-(a^3 ) Vo  p=l  &P  j£  ta-5^b21  c$y  ^s2,t2,vl^M2"Mlb21^ 

y  **o 


Hitter  Order  Approximation,  Higher  Dimensions. 


N 


.  n  a 

In  general,  if  EX  =  c  u  is  taken  as  the  basic  approximation. 


a=o 


la 


the  three-dimensional  integral  will  be  of  the  form 


2N  N  a 

<£>  s“  s5)  * 


where 


sa  3 


E  z 

a  W 


,(2) 


.(?) 

°ae> 


s~° '  =  E  E  E s- 


(3) 

'a&r 


,(P-D 


’  “a, a  ...a  , 

1  2  p-l 
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Thus,  increasing  accuracy  is  computationally  expensive,  and  increasing 
dimension  more  so,  but  nonetheless  the  process  should  compare  favorably 
with  numerical  integration,  as  borne  out  by  some  preliminary  comparisons, 
where  Romberg  integration  was  about  lOx  slower. 


11 


109 


2st^w^  (4tcw(5s2w+4  )-21 

945 


k=4 

Ut^w14  (5s2w+1 
945 


k=5 

Uat^w^ 

9^5 


k=0 

s  (t2  (w  (t2  (w  (t2  (b2\<  (2s2v (s2w+6  )+3  )-l4  )-21s2  (s2w+4  ) ) 

+  105)+315s2)-315)-315s2)/315 

k=l 

t2w(t2  (w  (t2  (s2w  (l0s2w  (s2w+4  )+l  )-12  )-63s2  (s2w+2  )  )+65 ) 

+  315s2)/315 


k^2 


lU 


-(8s8t6w5+56s6t6wl+-4sl4t6w5  -SUs6^*5  -192s2t6v2 


-  l»20s4t^v2^6t6w+lCK}8s2t^w+1260s4t2w-l89t4-l890a2t2 

-  9^5s4  )M5 

k=l 

Ust2w (t2 (w (t2 (s2w (2s2w (5s2w+24 )-J3  )-5U ) 

-  21s2 (3s2w+8 ) )+l89 )+315s2 )/915 


^t4w2(t2 (s2w(2082w(s2w+3 )-21)-9)-63s2(s2w+1 


-s (2s8t6w5  +l6s6t6w* -8akt6*?  -21s6t V -112s2t6w2 


il  ll  O  g  ,  O  h  h  O  _  h  „  Op 

-  126a  t  w  -*63t  w+462s  t  w+315s  t  w-l89t  -630a  t 

-  l89s1+)/l89 

k=l 

t2w(lOs8t4wl*+56s6t1+w5-36sl|t1+w2-63s6t2w2-l^s2t1+w 
-  210s4t2w+27tU+378s2t2 4515a1*  )/l89 


k=2 


_  st4w2 (2t2  (s2w  (2s2v ( 5s2w4-l8  )-23  )-l8  )-21s£: (3sS?+4 


2„_2 


159 


k=3 


W  (2t2(l0s2v(a2w4-2 


k=U 

2s^t8w4  (5s2y4-4 ) 


i  =  6 

k=0 

-  (4s10t^w'’+36s®t^wl4'  -36s^t^w^  -42s®t\^  -430s14 t^w2 

-  294s^tl|w2+450s2t^w+l470s!4t1,'w+630s^t2w-45t^-9U5a2tl+ 

-  1575s V -315s6 )/315 

k=l 

2st2w(l0s6t\^^4s6t\^-68sN;\*2-63  s^t2w2-300s2t1+w 
-  252s1;t2w+135t1++630s2t2+3l5s^)/315 

it  =5 

-  2s2t* w2 (2t2 (2s2w(s2w+5 ) (5s2w-4 )J*5 )-21s2 (3s2w+5 ) )/315 

k=3 

2s W  (2t2  (2s2w(5s2w+12  )-13  )-21s2 ) 

515 

k=4 

48^t^w*+(a2w+l) 

63 

k=5 

4s5t6v5 

'Ilf" 


J  =  7 

k  0 

a  (2s10t6w5+20s8t6w1|-31s6t6w5  -2is8tV  -366s4t6w2 
l68s6tl4w2+585s2t6w+1071s4tl4w+3l5s6t2w-135t6 
945s2t4-945«Ut2-135s6 )/l35 

k=l 

2t2w(1088t4w1>+72s6t1+w5-109s!+t^w2-63s6t2w2-51t0s2t1+w 

-  29Us4t2w+U05tU+945s2t2+315si+)/315 


k=2 


s^t^w2(t2(a2w(4a2w(3a2w^2U)-123  )-l80)»63s2(s2w+2)) 


k=3 

a^tV  (t2  (4s2w(3a2v-t-l^)-1*3  )-21a2 ) 
135 


k=3 


28^t*V4  (5s2w+6) 
135 


k=5 


2s6t6w 5 

tss” 
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lull 

k=0 

-  s2 (l6s10t6w5+l76s8t6w1*'-376s6t6w5 -l68s8tV 

-  4592s\6w2  -1512s6tlvw2+10080s2t6w+11760s1+t1+w 

+  2520s6t2w-3780t6-13230s2t*-8820s1*t2-945s6  )/945 

k=l 

8s5  t2w  (lOs8tUw1+  +80s6tV  -159s  \\»2  -68s6 t 2 w2 

-882s2tlvw-336s1;t2w+945t1++i323s2t2+315s1+)/945 

k=2 

Ss^t^w2  (t2  (s2w(4s2w(5s2w+27  )-175  )-315  )-2Is2(3s2w+7  )  )M5 

k=3 

8s^t\^  (t2(4s2w(5s2w+l6)-63)-21s2 
9^5 

k=4 

l6s6t6w^($B2w+7 ) 

9^5 

k=5 

l6s^t6w^ 

~~9%5 
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o  =  9 

k=0 

(2s10t^\?  +24s8t6w+  -^s^w5  -21s8tlvw5  -848sN;6w2 

210s6t\*2+2394s2t6w+1932s1+t\j+315s6t2v-126ot2 
26k6s2tl<'-1260s4t2-105s6  )/l05 


k=l 


4.2  8.4  4  -Q  6.4  3  „,o  4,4  2  6.2  2  2.4 

t  w(10s  t  w  +88s  t  w  -2l8s  t  w  -63s  t  w  -1344s  t  w 


-387s1+t2w+l890t4+lT64s2t2-t315s1<')/l05 


k=2 


A  -  10 


k=0 

-  s1+(4s10t6w5+52s8t6w1+-176s6t6w5-i+2s8tl4w5-239^sl+t6w2 

-  462s8t^w2+83l6s2t8w449l4s^t4’w+630s8t2w-5670t8 

-  7938s2t1;-2835s1+t2-l89s6)/i 89 


k=l 

2s^t2w  (10s8tl4’w)++96s8tl4’w^  -286sl4’ti+w2-63s8t2'w2 


-  1944s2t^w  -420s  **t2w+3  402 +2268s  2t2+3 15  s1*  )/l89 


k=2 

2s6t  V  (2t2  (s2w  (2s2w  (5s2w+33  )-155  )-378  )-63s2  (s2w+3 )  )/l89 


k=3 

2a^ (4t2 (5s2w(s2w+4)-27 )-21s2 ) 

k=4 

4s8t8w1+(5s2w+9) 

“  189 


k=5 

4.*tV 

”89 
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